This workflow detects the length of time for NDVI and EVI to go from baseline to peak over the course of the year. Each pixel is classified with a value reflecting the length of time in the year for NDVI and EVI to reach peak greenness.
Notes:
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3
library(rts)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(stringr)
library(foreach)
# Raster time series pulled from working directory + data/year/index
make_rts <- function(year = '2012', index = 'ndvi') {
# Just in case the year is passed as an integer
year <- as.character(year)
# **The data directory must be in the current working directory!
directory <- paste("data", year, index, "processed", sep='/')
# Get the GeoTIFF files
files <- list.files(directory, full.names=TRUE, pattern = "*.tif$")
# Pull all the files into a raster stack.
r <- stack(files)
# Extract dates from file names and cast as Date objects
dates <- as.Date(str_extract(files, '[0-9]{4}-[0-9]{2}-[0-9]{2}'))
# Reference the raster stack with their respective dates and collect them
# together into a time series
index_ts <- rts(r, dates)
# Give the raster time series result
return(index_ts)
}
days_to_max <- function(modists, ts_max, timepoints) {
start_of_year <- timepoints[1]
# Create an empty raster
reclass <- raster(crs = ts_max@crs, # Coordinate Ref Sys
ext = ts_max@extent, # Extent
resolution = res(ts_max), # Resolution
vals=NA) # Fill as NAs
# For each raster in TS...
foreach(i = 1:length(timepoints)) %do% {
# How many days in?
days_since <- as.numeric(timepoints[i] - start_of_year)
# Subtract the max from current raster in time series
test <- ts_max - modists[[i]]
# Redefine the 0 values as days since beginning of year
reclass[test == 0] <- days_since
}
return(reclass)
}
annual_stats_raster <- function(modists, FUN) {
# We have to use the nested index [[1]] to get the layer itself out of the rts object
return(apply.yearly(modists, FUN)[[1]])
}
# EVI and NDVI files for the year are moved in their own directory
# e.g. EVI files for 2012 are in data/2012/evi
# Get the time series
ndvi_ts <- make_rts(year = "2012", index = "ndvi")
evi_ts <- make_rts(year = "2012", index = "evi")
# Get rasters of statistics for each pixel over the course of the year
annual_max_ndvi <- annual_stats_raster(ndvi_ts, max)
annual_sd_ndvi <- annual_stats_raster(ndvi_ts, sd)
# Coefficient of variance calculation
# -- not working...?
#coeff_var <- function(x) sd(x)/mean(x)
#coeff_var_ndvi <- annual_stats_raster(ndvi_ts, coeff_var)
annual_max_evi <- annual_stats_raster(evi_ts, max)
annual_sd_evi <- annual_stats_raster(evi_ts, sd)
# These are the timepoints of our series
timepoints <- index(ndvi_ts@time)
reclass_ndvi <- days_to_max(ndvi_ts, annual_max_ndvi, timepoints)
reclass_evi <- days_to_max(evi_ts, annual_max_evi, timepoints)
# Dark violet to green
colors <- rev(rainbow(n = 300, s = 1, v = .7, start = .333, end = .8))
# Time series, multiplot
plot(ndvi_ts,
col=colors)
plot(evi_ts,
col=colors)
# Mean of all the rasters in each month,
plot(x = timepoints,
y = cellStats(ndvi_ts@raster, mean),
main = "2012 MODIS Raster Tiles\nMean NDVI Value Over Time",
type='o',
xlab = "month",
ylab = "mean NDVI")
plot(x = timepoints,
y = cellStats(evi_ts@raster, mean),
main = "2012 MODIS Raster Tiles\nMean EVI Value Over Time",
type='o',
xlab = "month",
ylab = "mean EVI")
plot(annual_sd_ndvi,
main = "2012 MODIS NDVI\nStandard Deviation for the Year",
col = colors)
plot(reclass_ndvi,
main = "2012 MODIS NDVI\nDays to Reach Maximum NDVI",
col = colors)
plot(annual_sd_evi,
main = "2012 MODIS EVI\nStandard Deviation for the Year",
col = colors)
plot(reclass_evi,
main = "2012 MODIS EVI\nDays to Reach Maximum EVI",
col = colors)